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ABSTRACT 

Conventional simulations of complex systems in the canonical ensemble suffer from 
the quasi-ergodicity problem. A simulation in generalized ensemble overcomes this diffi- 
culty by performing a random walk in potential energy space and other parameter space. 
From only one simulation run, one can obtain canonical-ensemble averages of physical 
quantities as functions of temperature by the single-histogram and/or multiple-histogram 
reweighting techniques. In this article we review the generalized-ensemble algorithms. 
Three well-known methods, namely, multicanonical algorithm, simulated tempering, and 
replica-exchange method, are described first. Both Monte Carlo and molecular dynamics 
versions of the algorithms are given. We then present further extensions of the above 
three methods. 



1 INTRODUCTION 

Canonical fixed-temperature simulations of complex systems such as spin glasses and 
biopolymers are greatly hampered by the multiple-minima problem, or the quasi-ergodicity 
problem. Because simulations at low temperatures tend to get trapped in a few of a huge 
number of local-minimum-energy states which are separated by high energy barriers, it is 
very difficult to obtain accurate canonical distributions at low temperatures by conven- 
tional Monte Carlo (MC) and molecular dynamics (MD) methods. One way to overcome 
this multiple-minima problem is to perform a simulation in a generalized ensemble where 
each state is weighted by an artificial, non-Boltzmann probability weight factor so that a 
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random walk in potential energy space may be realized (for reviews see, e.g., Refs. pQ 
[7]). The random walk allows the simulation to escape from any energy barrier and to 
sample much wider configurational space than by conventional methods. Monitoring the 
energy in a single simulation run, one can obtain not only the global-minimum-energy 
state but also canonical ensemble averages as functions of temperature by the single- 
histogram [HJ and/or multiple-histogram [91 [10] reweighting techniques (an extension of 
the multiple-histogram method is also referred to as weighted histogram analysis method 
(WHAM) [IQ]). Besides generalized-ensemble algorithms, which are usually based on local 
updates, methods based on non-local updates such as cluster algorithms and their gener- 
alizations have also been widely used [II] [13]. In this article, we focus our discussion on 
generalized-ensemble algorithms. 

One of the most well-known generalized-ensemble methods is perhaps multicanonical 
algorithm (MUCA) [HI EE] (for reviews see, e.g., Refs. [TIE 02]). (The method is also 
referred to as entropic sampling [18] and adaptive umbrella sampling [19J of the potential 
energy [20]. MUCA can also be considered as a sophisticated, ideal realization of a class 
of algorithms called umbrella sampling [21]. Also closely related methods are transition 
matrix methods reviewed in Refs. [221 S] an d random walk algorithm [231 121] > which is 
also referred to as density of states Monte Carlo [25] • See also Ref. [2E]-) MUCA and 
its generalizations have been applied to spin systems (see, e.g., Refs. [27J [32])- MUCA 
was also introduced to the molecular simulation field [33]. Since then MUCA and its 
generalizations have been extensively used in many applications in protein and related 
systems [34] [64]. Molecular dynamics version of MUCA has also been developed [4T1 1441 
[20] (see also Refs. [651 Sl| f° r Langevin dynamics version). MUCA has been extended so 
that flat distributions in other parameters instead of potential energy may be obtained 
[281 l29l l4"0"l |4"5"1 Wf\ [62] . Moreover, multidimensional (or multicomponent) extensions of 
MUCA can be found in Refs. [II S3 S3 EI] • 

While a simulation in multicanonical ensemble performs a free ID random walk in 
potential energy space, that in simulated tempering (ST) [6S1 E7J (the method is also 
referred to as the method of expanded ensemble |66j) performs a free random walk in 
temperature space (for a review, see, e.g., Ref. [68]). This random walk, in turn, induces 
a random walk in potential energy space and allows the simulation to escape from states of 
energy local minima. ST has also been applied to protein folding problem [691 S21 S21 [70]. 

The generalized-ensemble algorithms are powerful, but in the above two methods 
the probability weight factors are not a priori known and have to be determined by 
iterations of short trial simulations. This process can be non-trivial and very tedius for 
complex systems with many degreees of freedom. Therefore, there have been attempts to 
accelerate the convergence of the iterative process for MUCA weight factor determination 
[23 SQl EU [721 (731 120] (see also Refs. [161 El]). 

In the replica- exchange method (REM) [75] [77], the difficulty of weight factor deter- 
mination is greatly alleviated. (A closely related method was independently developed in 
Ref. [78J. Similar methods in which the same equations are used but emphasis is laid on 
optimizations have been developed [791 EH] . REM is also referred to as multiple Markov 
chain method [81] and parallel tempering |68j. Details of literature about REM and re- 
lated algorithms can be found in recent reviews [521 13 ■) I n this method, a number of 
non-interacting copies (or replicas) of the original system at different temperatures are 
simulated independently and simultaneously by the conventional MC or MD method. Ev- 
ery few steps, pairs of replicas are exchanged with a specified transition probability. The 
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weight factor is just the product of Boltzmann factors, and so it is essentially known. 

REM has already been used in many applications in protein systems [531 EH 170] [85]- 
|97j . Other molecular simulation fields have also been studied by this method in various 
ensembles [98j [103] . Moreover, REM was applied to cluster studies in quantum chemistry 
field [104] . The details of molecular dynamics algorithm have been worked out for REM 
in Ref. [M] (see also Refs. [831 HOlj ) ■ This led to a wide application of replica-exchange 
molecular dynamics method in the protein folding problem |105j - [TT2] . 

However, REM also has a computational difficulty: As the number of degrees of 
freedom of the system increases, the required number of replicas also greatly increases, 
whereas only a single replica is simulated in MUCA or ST. This demands a lot of computer 
power for complex systems. Our solution to this problem is: Use REM for the weight fac- 
tor determinations of MUCA or ST, which is much simpler than previous iterative methods 
of weight determinations, and then perform a long MUCA or ST production run. The first 
example is the replica- exchange multicanonical algorithm (REMUCA) [55 1 1931 IM] . In RE- 
MUCA, a short replica-exchange simulation is performed, and the multicanonical weight 
factor is determined by the multiple-histogram reweighting techniques [HI [10]. Another 
example of such a combination is the replica- exchange simulated tempering (REST) [89] . 
In REST, a short replica-exchange simulation is performed, and the simulated tempering 
weight factor is determined by the multiple-histogram reweighting techniques [HI dO] • 

We have introduced two further extensions of REM, which we refer to as multicanon- 
ical replica- exchange method (MUCAREM) [551 ESI EI] (see also Refs. [1131 HI]) and 
simulated tempering replica- exchange method (STREM) [115] (see also Ref. [116] for a 
similar idea). In MUCAREM, a replica-exchange simulation is performed with a small 
number of replicas each in multicanonical ensemble of different energy ranges. In STREM, 
on the other hand, a replica-exchange simulation is performed with a small number of 
replicas in "simulated tempering" ensemble of different temperature ranges. 

Finally, one is naturally led to a multidimensional (or, multivariable) extension of 
REM, which we refer to as multidimensional replica- exhcange method (MREM) [86] (see 
also Refs. [TTT1 EH1 ESI E21 EH]). A special realization of MREM is replica- exchange 
umbrella sampling (REUS) [86] and it is particularly useful in free energy calculations (see 
also Ref. [57] for a similar idea). 

In this article, we describe the generalized-ensemble algorithms mentioned above. 
Namely, we first review the three familiar methods: MUCA, ST, and REM. We then 
present further extensions of the three methods. 

2 GENERALIZED-ENSEMBLE ALGORITHMS 

2.1 Multicanonical Algorithm and Simulated Tempering 

Let us consider a system of N atoms of mass (k = 1, • • • , N) with their coordinate 
vectors and momentum vectors denoted by q = {qfi, • • • , qfjv} an d p = {Pi, • • • ,Pn}> 
respectively. The Hamiltonian H(q,p) of the system is the sum of the kinetic energy 
K{p) and the potential energy E(q): 

H(q,p) = K(p) + E(q) , (1) 



3 



where 

N 2 

Kip) = T,P~- (2) 

In the canonical ensemble at temperature T each state x = (q,p) with the Hamiltonian 
H(q,p) is weighted by the Boltzmann factor: 

W B (x;T) = exp(-l3H(q,p)) , (3) 

where the inverse temperature (3 is defined by (3 = l/k B T (k B is the Boltzmann constant). 
The average kinetic energy at temperature T is then given by 

(^)>t=(e^ =!W. (4) 



^ 2m k I 2 



Because the coordinates q and momenta p are decoupled in Eq. (CQ), we can suppress 
the kinetic energy part and can write the Boltzmann factor as 

W B (x; T) = W B (E; T) = exp(-/3E) . (5) 

The canonical probability distribution of potential energy P B (E;T) is then given by the 
product of the density of states n(E) and the Boltzmann weight factor W B (E; T): 

P B (E;T)<xn(E)W B (E;T) . (6) 

Since n{E) is a rapidly increasing function and the Boltzmann factor decreases expo- 
nentially, the canonical ensemble yields a bell-shaped distribution which has a maximum 
around the average energy at temperature T. The conventional MC or MD simulations 
at constant temperature are expected to yield P B (E; T). A MC simulation based on the 
Metropolis algorithm [120] is performed with the following transition probability from a 
state x of potential energy E to a state x' of potential energy E'\ 

( W B (E'-T)\ 

w(x -> x ) = min I 1, ^— ^i— 1 = min (1, exp (—/3AE)) . (7) 

where 

AE — E' — E . (8) 

A MD simulation, on the other hand, is based on the following Newton equations of 
motion: 

q k = — , (9) 

m k 

dE 

P k = - Wh =fk> ( 10 ) 

where f k is the force acting on the k-th atom [k = 1, ■ ■ • , N). This set of equations actually 
yield the microcanonical ensemble, and we have to add a thermostat in order to obtain the 
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canonical ensemble at temperature T. Here, we just follow Nose's prescription |121l 1122] . 
and we have 

q k = — , (ii) 

dE s s 

p - = -W k '~s Pk = f *~i p - (12 > 

* = I") 

N 2 

Ps = J2^-3Nk B T = 3Nk B (T(t)-T) , (14) 



k=i 



where s is Nose's scaling parameter, Q is its mass, P s is its conjugate momentum, and 
the "instantaneous temperature" T(t) is defined by 

However, in practice, it is very difficult to obtain accurate canonical distributions of 
complex systems at low temperatures by conventional MC or MD simulation methods. 
This is because simulations at low temperatures tend to get trapped in one or a few of 
local-minimum-energy states. 

In the multicanonical ensemble [TJ1 [15], on the other hand, each state is weighted by 
a non-Boltzmann weight factor W mu (E) (which we refer to as the multicanonical weight 
factor) so that a uniform potential energy distribution P mn (E) is obtained: 

P mu (E) oc n(E)W mu (E) = const . (16) 

The flat distribution implies that a free random walk in the potential energy space is real- 
ized in this ensemble. This allows the simulation to escape from any local minimum-energy 
states and to sample the configurational space much more widely than the conventional 
canonical MC or MD methods. 

The definition in Eq. f[T6|) implies that the multicanonical weight factor is inversely 
proportional to the density of states, and we can write it as follows: 

W mn (E) = exp [~A)£mu(£; To)] = , (17) 

where we have chosen an arbitrary reference temperature, T = l/k^fio, and the "multi- 
canonical 'potential energy" is defined by 

E mvL (E; T ) = k B T Q \nn(E) = T S(E) . (18) 

Here, S(E) is the entropy in the microcanonical ensemble. Since the density of states of 
the system is usually unknown, the multicanonical weight factor has to be determined 
numerically by iterations of short preliminary runs [HI US] . 

A multicanonical MC simulation is performed, for instance, with the usual Metropolis 
criterion |120] : The transition probability of state x with potential energy E to state x' 
with potential energy E' is given by 

w(x -> x') = min (l, ) = min (l, = min (1, exp {-fa AE mu )) , (19) 
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where 

AE mu = E ma (E'; T ) - E mu (E; T ) . (20) 

The MD algorithm in the multicanonical ensemble also naturally follows from Eq. (TTTj) . 
in which the regular constant temperature MD simulation (with T = T ) is performed by 
replacing E by E mn in Eq. (JT2D [III 



dE mu (E;T ) s dE mu (E; T ) a 

Pfc = ^ S Pk = dE h-,*' (21) 

From Eq. (ITS]) this equation can be rewritten as 

Pk = f k -~ s Pk- ( 22 ) 

where the following thermodynamic relation gives the definition of the "effective temper- 
ature" T(E): 

9 -^l - -J- (23) 

with 

E a = <E> T{Ea) . (24) 

If the exact multicanonical weight factor W mu (E) is known, one can calculate the 
ensemble averages of any physical quantity A at any temperature T (= as follows: 

£ A(£)P B (£;T) £ A(£7)n(£7)exp(-^) 

< A > r = = -L , (25) 

^ Pb(^;T) £ n{E)e*p{-ffl 

E E 

where the density of states is given by (see Eq. (TT?j) ) 

= • (26) 

The summation instead of integration is used in Eq. (125]) . because we often discretize the 
potential energy E with step size e (i? = i = 1, 2, ■ ■ •). Here, the explicit form of the 
physical quantity A should be known as a function of potential energy E. For instance, 
A(E) = E gives the average potential energy < E >t as a function of temperature, and 
A(E) = (3 2 (E- < E > T ) 2 gives specific heat. 

In general, the multicanonical weight factor W mn (E), or the density of states n(E), is 
not a priori known, and one needs its estimator for a numerical simulation. This estimator 
is usually obtained from iterations of short trial multicanonical simulations. The details 
of this process are described, for instance, in Refs. [271 [36] • However, the iterative process 
can be non-trivial and very tedius for complex systems. 

In practice, it is impossible to obtain the ideal multicanonical weight factor with com- 
pletely uniform potential energy distribution. The question is when to stop the iteration 
for the weight factor determination. Our criterion for a satisfactory weight factor is that 
as long as we do get a random walk in potential energy space, the probability distribution 
Pmu{E) does not have to be completely flat with a tolerance of, say, an order of magnitude 
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deviation. In such a case, we usually perform with this weight factor a multicanonical 
simulation with high statistics (production run) in order to get even better estimate of the 
density of states. Let N mu (E) be the histogram of potential energy distribution P mu (E) 
obtained by this production run. The best estimate of the density of states can then be 
given by the single-histogram reweighting techniques [8] as follows (see the proportionality 
relation in Eq. ffTBl ): 

= t§| ' (27) 

By substituting this quantity into Eq. (1251) . one can calculate ensemble averages of phys- 
ical quantity A(E) as a function of temperature. Moreover, ensemble averages of any 
physical quantity A (including those that cannot be expressed as functions of potential 
energy) at any temperature T (= \/k-Q0) can now be obtained as long as one stores the 
"trajectory" of configurations (and A) from the production run. Namely, we have 

no 



J2A(x(k))W^(E(x(k)))exp[-PE(x(k))} 

no 

J2W-^(E(x(k)))exp[-f3E(x(k))] 



< A > T = k -^r n , (28) 



k=l 



where x(k) is the configuration at the k-th MC (or MD) step and Uq is the total number of 
configurations stored. Note that when A is a function of E, Eq. (1281) reduces to Eq. (1251) 
where the density of states is given by Eq. (1271) . 

Eqs. ( 1251) and ( 1281) or any other equations which involve summations of exponential 
functions often encounter with numerical difficulties such as overflows. These can be 
overcome by using, for instance, the following equation [1231 H24j : For C = A + B (with 
A > and B > 0) we have 



InC =ln 



max(A, B) 



min(A, B) 
max(A, B) 



(29) 



max(lnyl, In I?) + In {1 + exp [min(lnA, In B) — max (In A, In B)]} 



We now briefly review the original simulated tempering (ST) method [66l [67] . In this 
method temperature itself becomes a dynamical variable, and both the configuration and 
the temperature are updated during the simulation with a weight: 

W S t(E; T) = exp {-(3E + a{T)) , (30) 

where the function a(T) is chosen so that the probability distribution of temperature is 
flat: 

P ST (T) =JdE n(E) W S t(E; T) = J dE n(E) exp {—f3E + a(T)) = const . (31) 

Hence, in simulated tempering the temperature is sampled uniformly. A free random walk 
in temperature space is realized, which in turn induces a random walk in potential energy 
space and allows the simulation to escape from states of energy local minima. 

In the numerical work we discretize the temperature in M different values, T m (m = 
1, ■ • • , M). Without loss of generality we can order the temperature so that T\ < T2 < 
■ ■ ■ < Tm. The lowest temperature T\ should be sufficiently low so that the simulation 



7 



can explore the global-minimum-energy region, and the highest temperature Tm should 
be sufficiently high so that no trapping in an energy-local-minimum state occurs. The 
probability weight factor in Eq. (1301) is now written as 

W S t{E; T m ) = exp{-(3 m E + a m ) , (32) 

where a m = a(T m ) (m — 1, ■ • • , M). Note that from Eqs. fl3~Tl) and fl32l we have 

exp(— a m ) oc J dE n(E) exp(—f3 rn E) . (33) 

The parameters a m are therefore "dimensionless" Helmholtz free energy at temperature 
T m (i.e., the inverse temperature (3 m multiplied by the Helmholtz free energy). We re- 
mark that the density of states n(E) (and hence, the multicanonical weight factor) and 
the simulated tempering weight factor a m are related by a Laplace transform [12]. The 
knowledge of one implies that of the other, although in numerical work the inverse Laplace 
transform of Eq. (133]) is nontrivial. 

Once the parameters Clfn £1X6 determined and the initial configuration and the initial 
temperature T m are chosen, a simulated tempering simulation is then realized by alter- 
nately performing the following two steps [HE1 EZ] : 

1. A canonical MC or MD simulation at the fixed temperature T m (based on Eq. (|7j) 
or Eq. (flOj) ) is carried out for a certain steps. 

2. The temperature T m is updated to the neighboring values T m±1 with the configura- 
tion fixed. The transition probability of this temperature-updating process is given 
by the Metropolis criterion (see Eq. (132]) ): 



w(T m -> T m±1 ) = min |l, ^^-^y) = min ( X > ex P (~ A )) > ( 34 ) 

where 

A = (/3 m ±i - f3 m ) E - (a m ±i - a m ) . (35) 

Note that in Step 2 we exchange only pairs of neighboring temperatures in order to secure 
sufficiently large acceptance ratio of temperature updates. 

As in multicanonical algorithm, the simulated tempering parameters a m = a(T m ) 
(m = 1, ■ ■ • , M) are also determined by iterations of short trial simulations (see, e.g., 
Refs. [681 [69], H3] f° r details). This process can be non-trivial and very tedius for complex 
systems. 

After the optimal simulated tempering weight factor is determined, one performs a 
long simulated tempering run once. The canonical expectation value of a physical quantity 
A at temperature T m (m = 1, • • • , M) can be calculated by the usual arithmetic mean as 
follows: 

n m 

<A> Tm =—J2A(x m (k)) , (36) 

where x m (k) (k — 1, • ■ • ,n m ) are the configurations obtained at temperature T m and n m 
is the total number of measurements made at T = T m . The expectation value at any 
intermediate temperature can also be obtained from Eq. (125]) . where the density of states 
is given by the multiple-histogram reweighting techniques [9] [10] as follows. Let N m (E) 
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and n m be respectively the potential-energy histogram and the total number of samples 
obtained at temperature T m = l/ksPm (m = 1, • • • , M). The best estimate of the density 
of states is then given by [HI CD] 

M 

E 9 m l N m (E) 

n{E) = - m=1 , (37) 

E 9m n rn exp(/ m - (3 m E) 
m=l 

where we have for each m (= 1, • ■ • , M) 

exp(-/ m ) = E n(E) exp(-fJ m E) . (38) 

E 

Here, g m = 1 + 2r m , and r m is the integrated autocorrelation time at temperature T m . 
For many systems the quantity g m can safely be set to be a constant in the reweighting 
formulae [10], and hereafter we set g m — 1. 

Note that Eqs. ( !37l) and ( |38l) are solved self-consistently by iteration [9], [10] to obtain 
the density of states n(E) and the dimensionless Helmholtz free energy f m . Namely, we 
can set all the f m (m = 1, • • • , M) to, e.g., zero initially. We then use Eq. ( 1371) to obtain 
n(E), which is substituted into Eq. (|38|) to obtain next values of f m , and so on. 

Moreover, ensemble averages of any physical quantity A (including those that cannot 
be expressed as functions of potential energy) at any temperature T (= l/fce/?) can now 
be obtained from the "trajectory" of configurations of the production run. Namely, we 
first obtain f m (m = 1, • • • , M) by solving Eqs. (13T1) and (13"8"|) self-consistently, and then 
we have [93] 

M n m i 

£ 5>M*0) m ex P [-PE(x m (k))} 

J2n e exp lf e - (3 e E(x m (k))} 

< A >t= inz 1 ' ( 39 ) 

E E ~m ex P [-PE{x m {k))\ 

m-l fc=i j2 ne exp [f g _ /3^( Xm (A;))] 
where x m {k) {k = 1, • ■ ■ ,n m ) are the configurations obtained at temperature T m . 



2.2 Replica-Exchange Method 

The replica- exchange method (REM) [75]-[77] was developed as an extension of simulated 
tempering [75] (thus it is also referred to as parallel tempering [68]) (see, e.g., Ref. [HI] 
for a detailed description of the algorithm). The system for REM consists of M non- 
interacting copies (or, replicas) of the original system in the canonical ensemble at M 
different temperatures T rn (m — 1, ■ • • , M). We arrange the replicas so that there is always 
exactly one replica at each temperature. Then there exists a one-to-one correspondence 
between replicas and temperatures; the label % {% = 1, • • • , M) for replicas is a permutation 
of the label m (m = 1, • • • , M) for temperatures, and vice versa: 

| i = i( m ) = f{m) , , 40 > 
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where f(m) is a permutation function of m and f~ l {i) is its inverse. 

Let X = , ■ • • ,x{^ M ^| = j^^!), • • ■ ;^m(if)} stand for a "state" in this general- 

ized ensemble. Each "substate" is specified by the coordinates and momenta p^ 
of N atoms in replica i at temperature T m : 

xg=( g M,pM) . (41) 



Because the replicas are non-interacting, the weight factor for the state X in this 
generalized ensemble is given by the product of Boltzmann factors for each replica (or at 
each temperature): 



WremPO =f[exp{-(3 m{i) H (qW,pW)} = f[ exp {-/3 m H (V^/Ml 



i=l m=l 

{M 1 r M 

-J2Pm(i)H (V* 1 ,^ 1 ) = exp - J2 PmH (q 
i=l ) I m=l 



"Ml ^M" 1 )] 



(42) 



where i(m) and m(z) are the permutation functions in Eq. (1401) . 

We now consider exchanging a pair of replicas in the generalized ensemble. Suppose 
we exchange replicas i and j which are at temperatures T m and T n , respectively: 

X — |- • • , xQ, ■ ■ • , x]j\ - • •} > X — |- • ■ , x]jQ , ■ ■ • , , • ■ -| . (43) 

Here, i, j, m, and n are related by the permutation functions in Eq. (1401) . and the exchange 
of replicas introduces a new permutation function /': 

% = f(m) — > j = f'(m) , 
{ ■ , \ f ,\ 44 

1 J = f{n) — ► % = fin) . 



The exchange of replicas can be written in more detail as 

') . 

(45) 



where the definitions for p^' and p 1 ^' will be given below. We remark that this process is 
equivalent to exchanging a pair of temperatures T m and T n for the corresponding replicas 
i and j as follows: 

T [i] = f/tH n [i]\ > r [i]' = f/yW r)W /N ) 

In the original implementation of the replica- exchange method (REM) [75] [77], Monte 
Carlo algorithm was used, and only the coordinates q (and the potential energy function 
E(q)) had to be taken into account. In molecular dynamics algorithm, on the other 
hand, we also have to deal with the momenta p. We proposed the following momentum 
assignment in Eq. (|4"5]1 (and in Eq. (T461) ) [84] : 

(47) 



pW = \l — — p 





p ui =\i — p 
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which we believe is the simplest and the most natural. This assignment means that we 
just rescale uniformly the velocities of all the atoms in the replicas by the square root of 
the ratio of the two temperatures so that the temperature condition in Eq. (J4J) may be 
satisfied. 

In order for this exchange process to converge towards an equilibrium distribution, 
it is sufficient to impose the detailed balance condition on the transition probability 
w{X^X')\ 



w REM (x) 



w(X -> X') 



Wrem(X') 



w(X' -> X) 



(48) 



Z y ' Z 

where Z is the partition function of the entire system. From Eqs. ((H), (j2J), ( )42l) . (J47j), and 
(1481) . we have 

W REM (X) 



exp {-/3 m [K (pW) + E (qW)] - n [K (pW) + £ (gW) 

+/3 m jK (pW) + £ (gl^] + /?„ [# (p b1 ) + E (q®)] } 
exp {-/3 m7 ^X (p bl ) - A^if (pW) + /3 m # (pW) + /3 n # 
'E (qW) - E (gW)] - /3„ (gW) - £ (g^)] } , 



P 



li] 



(49) 



exp (-A) 



where 



E(q 



Pm (E (q 



Pn(E(q 



E (gW) 



(50) 
(51) 



and i, j, m, and n are related by the permutation functions in Eq. (140]) before the exchange: 

(52) 



This can be satisfied, for instance, by the usual Metropolis criterion |120j (see also Eqs. ([7]), 
(USD, and (|MD): 



w(X -> X') 



W X 



X 



[j] 



min (1, exp (—A)) 



(53) 



where in the second expression (i.e., u>(a;Jj}|a;^)) we explicitly wrote the pair of replicas 
(and temperatures) to be exchanged. Note that this is exactly the same criterion that 
was originally derived for Monte Carlo algorithm [75]-|77j. 

Without loss of generality we can again assume T\ < T 2 < ■ ■ ■ < Tm- A simulation of 
the replica- exchange method (REM) [7S] [77] is then realized by alternately performing 
the following two steps: 

1. Each replica in canonical ensemble of the fixed temperature is simulated simultaneously 
and independently for a certain MC or MD steps. 

2. A pair of replicas at neighboring temperatures, say xfy and x^} +1 , are exchanged 
with the probability w (x^ %m+i) in Eq. (|53|) . 
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Note that in Step 2 we exchange only pairs of replicas corresponding to neighboring tem- 
peratures, because the acceptance ratio of the exchange process decreases exponentially 
with the difference of the two /3's (see Eqs. ( I5"lj) and (1351) ). Note also that whenever a 
replica exchange is accepted in Step 2, the permutation functions in Eq. (j4"Ul) are updated. 

The REM simulation is particularly suitable for parallel computers. Because one can 
minimize the amount of information exchanged among nodes, it is best to assign each 
replica to each node (exchanging pairs of temperature values among nodes is much faster 
than exchanging coordinates and momenta). This means that we keep track of the per- 
mutation function m(i; t) = / _1 (i; t) in Eq. (1401) as a function of MC or MD step t during 
the simulation. After parallel canonical MC or MD simulations for a certain steps (Step 
1), M/2 pairs of replicas corresponding to neighboring temperatures are simulateneously 
exchanged (Step 2), and the pairing is alternated between the two possible choices, i.e., 
(T\,T 2 ), (T 3 ,T 4 ), ■•• and (T 2 ,T 3 ), (T 4 ,T 5 ), ••, 

The major advantage of REM over other generalized-ensemble methods such as mul- 
ticanonical algorithm p3J [15] and simulated tempering [661 EZ] lies in the fact that the 
weight factor is a priori known (see Eq. (]42p ). while in the latter algorithms the deter- 
mination of the weight factors can be very tedius and time-consuming. A random walk 
in "temperature space" is realized for each replica, which in turn induces a random walk 
in potential energy space. This alleviates the problem of getting trapped in states of 
energy local minima. In REM, however, the number of required replicas increases as the 
system size N increases (according to \^N) [75]. This demands a lot of computer power 
for complex systems. 

2.3 Replica-Exchange Multicanonical Algorithm and Replica- 
Exchange Simulated Tempering 

The replica- exchange multicanonical algorithm (REMUCA) [881 [931 E3] overcomes both 
the difficulties of MUCA (the multicanonical weight factor determination is non-trivial) 
and REM (a lot of replicas, or computation time, is required). In REMUCA we first 
perform a short REM simulation (with M replicas) to determine the multicanonical weight 
factor and then perform with this weight factor a regular multicanonical simulation with 
high statistics. The first step is accomplished by the multiple-histogram reweighting 
techniques [10] - Let N m (E) and n m be respectively the potential-energy histogram and 
the total number of samples obtained at temperature T m (= I/ZceA™) of the REM run. 
The density of states n(E) is then given by solving Eqs. (1571) and (1381) self-consistently by 
iteration. 

Once the estimate of the density of states is obtained, the multicanonical weight 
factor can be directly determined from Eq. (ITT]) (see also Eq. ( ITS]) ). Actually, the density 
of states n(E) and the multicanonical potential energy, E mu (E;To), thus determined are 
only reliable in the following range: 

E x <E<E Ml (54) 

where 




<E> n , 
<E>t m , 



(55) 
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and Ti and T M are respectively the lowest and the highest temperatures used in the REM 
run. Outside this range we extrapolate the multicanonical potential energy linearly: [HH] 




(E — Ei) + E mu (Ei, T ) 



E=E X 



(E — E M ) + E mu (E M ; T Q 



for E < Ei, 

for Ei<E< E M , 

for E > E M . 



(56) 



E=E M 



The multicanonical MC and MD runs are then performed respectively with the Metropolis 
criterion of Eq. (fT9|) and with the modified Newton equation in Eq. (12~T1) . in which £^(E) 
in Eq. (156]) is substituted into E mu (E;T ). We expect to obtain a flat potential energy 
distribution in the range of Eq. (1541) . Finally, the results are analyzed by the single- 
histogram reweighting techniques as described in Eq. (1271) (and Eq. ([25]) ). 

Some remarks are now in order. From Eqs. (1181) . (B3"j) . (1241) , and (I55p . Eq. (1561) becomes 



(E-E 1 )+T S(Ei 



(To 

Ti 

T S(E) , 
^{E - E M 

J-M 



T, 



Ti 



E + const 



+ ToS(E M ) = -^E + const , for E > E M 

J-M 



for E < Ei =< E > Tl , 
for Ei<E< Em, 

< E >T M - 



The Newton equation in Eq. ([2TD is then written as (see Eqs. (I22j) . (1231). and fl24j) ; 

for E < Ei, 



(57) 



Pk 



To 

rp J k 

¥ 

-t 

T(E) 
% 

Ti J A: 
M 



f 



s 

- Pk 

s 

s 



Pk 



s 

- Pk 

s 



for Ei<E< E M , 
for E > E M - 



(58 



Because only the product of inverse temperature /3 and potential energy E enters in the 
Boltzmann factor (see Eq. a rescaling of the potential energy (or force) by a constant, 
say a, can be considered as the rescaling of the temperature by 1/a [4T| 1101] . Hence, our 
choice of S^(E) in Eq. ([56]) results in a canonical simulation at T — Ti for E < Ei, a 
multicanonical simulation for E\ < E < i?jv/ ; and a canonical simulation at T = for 
> -Em. Note also that the above arguments are independent of the value of T , and we 
will get the same results, regardless of its value. 

For Monte Carlo method, the above statement follows directly from the following 
equation. Namely, our choice of the multicanonical potential energy in Eq. (1561) gives (by 
substituting Eq. (J57J into Eq. 071) ) 



W mu (E) =exp\-p £W(E) 



exp {—(3iE + const) 
1 



n(E) ' 

exp {—(3 M E + const) 



for E < Ei, 

for Ei<E< E M , 

for E > E M . 



(59) 



We now present another effective method of the multicanonical weight factor determi- 
nation [3], which is closely related to REMUCA. We first perform a short REM simula- 
tion as in REMUCA and calculate < E >t as a function of T by the multiple-histogram 
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reweighting techniques (see Eqs. (157]) and ([38]) ). Let us recall the Newton equation of 
Eq. (1221) and the thermodynamic relation of Eqs. (1231) and ([M]). The effective temperature 



REM simulation by the multiple-histogram reweighting techniques. Given its derivative, 
the multicanonical potential energy can then be obtained by numerical integration (see 



We remark that the same equation was used to obtain the multicanonical weight factor in 
Ref. [72], where < E > T was estimated by simulated annealing instead of REM. Essen- 
tially the same formulation was also recently used in Ref. [61] to obtain the multicanonical 
potential energy, where < E was calculated by conventional canonical simulations. 

We finally present the new method which we refer to as the replica- exchange simulated 
tempering (REST) [SH]. In this method, just as in REMUCA, we first perform a short 
REM simulation (with M replicas) to determine the simulated tempering weight factor 
and then perform with this weight factor a regular ST simulation with high statistics. 
The first step is accomplished by the multiple-histogram reweighting techniques [9l HP], 
which give the dimensionless Helmholtz free energy f m (see Eqs. ( 1371) and ( [38]) ). 

Once the estimate of the dimensionless Helmholtz free energy f m are obtained, the 
simulated tempering weight factor can be directly determined by using Eq. (132]) where we 
set a m = f m (compare Eq. ([33]) with Eq. ([38]) ). A long simulated tempering run is then 
performed with this weight factor. Let N m (E) and n m be respectively the potential-energy 
histogram and the total number of samples obtained at temperature T m (= l/k B /3 m ) from 
this simulated tempering run. The multiple-histogram reweighting techniques of Eqs. (I3"7]) 
and (138]) can be used again to obtain the best estimate of the density of states n(E). The 
expectation value of a physical quantity A at any temperature T (= l/ks/3) is then 
calculated from Eq. ([25]) . 

The formulations of REMUCA and REST are simple and straightforward, but the 
numerical improvement is great, because the weight factor determination for MUCA and 
ST becomes very difficult by the usual iterative processes for complex systems. 

2.4 Multicanonical Replica- Exchange Method and Simulated Tem- 
pering Replica-Exchange Method 

In the previous subsection we presented REMUCA, which uses a short REM run for the 
determination of the multicanonical weight factor. Here, we present two modifications 
of REM and refer the new methods as multicanonical replica-exchange method (MU- 
CAREM) [SSI ESI El] and simulated tempering replica-exchange method (STREM) [TIB] . 
In MUCAREM the production run is a REM simulation with a few replicas not in the 
canonical ensemble but in the multicanonical ensemble, i.e., different replicas perform 
MUCA simulations with different energy ranges. Likewise in STREM the production run 
is a REM simulation with a few replicas that performs ST simulations with different tem- 
perature ranges. While MUCA and ST simulations are usually based on local updates, 
a replica-exchange process can be considered to be a global update, and global updates 
enhance the sampling further. 




Eqs. (US]) and ([23])): [3] 




(60) 
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We first describe MUCAREM. Let A4 be the number of replicas. Here, each replica is 
in one-to-one correspondence not with temperature but with multicanonical weight factors 
of different energy range. Note that because multicanonical simulations cover much wider 
energy ranges than regular canonical simulations, the number of required replicas for the 
production run of MUCAREM is much less than that for the regular REM (M. <C M). 
The weight factor for this generalized ensemble is now given by (see Eq. ()42p ) 

M M 
i=l m=l 



WWarem(A) = J] Wi™ »} (E = J] ^' (E (xg™)])) , (61) 



where we prepare the multicanonical weight factor (and the density of states) separately 
for m regions (see Eq. (fT7])): 

(E (sg)) = exp [-(3 m £™ (E (xg))] = ^ * fJ ^ . (62) 



Here, we have introduced M. arbitrary reference temperatures T m = l/k^f3 m (m = 
but the final results will be independent of the values of T m , as one can 
see from the second equality in Eq. (162~1) (these arbitrary temperatures are necessary only 
for MD simulations). 

Each multicanonical weight factor W^'(E), or the density of states n^ m '(E), is 
defined as follows. For each m (m = 1, ■■■,A4), we assign a pair of temperatures 
(r L W ,rf). Here, we assume that T L {m} < Tt } and arrange the temperatures so 
that the neighboring regions covered by the pairs have sufficient overlaps. Without loss 
of generality we can assume < ■ ■ ■ < T^ M ^ and < ■ ■ ■ < . We define the 
following quantities: 

Et } = <E> TLim} , 
EjT } = <E> Tn{m} , (m = l,---,M). 

Suppose that the multicanonical weight factor W mu (E) (or equivalently, the multi- 
canonical potential energy E mu (E;T Q ) in Eq. (fTBj) ) has been obtained as in REMUCA or 
by any other methods in the entire energy range of interest (E^ < E < Er 4 ). We then 
have for each m (m = 1, • • • , Ai) the following multicanonical potential energies (see Eq. 

my- m 



dE mn (E^ '; T, 



dE 



m) (E - Et } ) + E mu (E[ m} ; T m ) , for £ < Et } 



E mn (E; T m ) , for Et } <E< E { K m \ (64) 



<9£ 



£ - E^ + £ mu *T ' ; T m , for £ > E£ 



Finally, a MUCAREM simulation is realized by alternately performing the following 
two steps. 

1. Each replica of the fixed multicanonical ensemble is simulated simultaneously and 
independently for a certain MC or MD steps. 
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2. A pair of replicas, say i and j, which are in neighboring multicanonical ensembles, 
say m-th and (m+l)-th, respectively, are exchanged: X = {• • • , xM, • • • , x^ +1 , ■ ■ •} — 

X' = {• • • , x^j, • • • , • • •}. The transition probability of this replica exchange 
is given by the Metropolis criterion: 

w(X^X') =min(l,exp(-A)) , (65) 

where we now have (see Eq. (1501) ) [88] 



A =A-{£fe } (^(ff w ))-* } (^(? w ))} 

- /w {4r 1} - 4r 1} (^ ] ))} 



(66) 



Here, E (a^) and .E (q^ are the potential energy of the i-th replica and the j-th 
replica, respectively 

Note that in Eq. ( 1661) we need to newly evaluate the multicanonical potential energy, 
£^ } (E(q^)) and £^ +1} (E(q^)), because £^ } (E) and £^( E ) are, in general, different 
functions for m ^ n. 

In this algorithm, the m-th multicanonical ensemble actually results in a canonical 
simulation at T = T^ m ^ for E < E^, a multicanonical simulation for E^ < E < E^, 
and a canonical simulation at T = T^ 71 ^ for E > E^™\ while the replica-exchange process 
samples states of the whole energy range (E^ < E < E^^). 

For obtaining the canonical distributions at any intermediate temperature T, the 
multiple-histogram reweighting techniques [H [10] are again used. Let N m (E) and n m 
be respectively the potential-energy histogram and the total number of samples obtained 
with the multicanonical weight factor W^{E) (m = 1, • ■ ■ , A4). The expectation value 
of a physical quantity A at any temperature T (= 1//cb/3) is then obtained from Eq. ( 1251) . 
where the best estimate of the density of states is obtained by solving the WHAM equa- 
tions, which now read [88] 

M M 

N m (E) N m(E) 

n(E) = M m=1 = ^ ^ , (67) 

n m eMfm)WW(E) J2 n m exp (f m - (3 m £^{E) 

m=l m=l 

and for each m (= 1, • ■ ■ , M) 

eM-U) = Y<E)W^\E) = £ n(E) exp (-P m £^ } (E)) . (68) 

E E 

Note that W^(E) is used instead of the Boltzmann factor exp(— (3 m E) in Eqs. (1371) and 

(EHD- 

Moreover, ensemble averages of any physical quantity A (including those that cannot 
be expressed as functions of potential energy) at any temperature T (= l/k^fi) can now 
be obtained from the "trajectory" of configurations of the production run. Namely, we 
first obtain f m (m — 1, • • • ,M) by solving Eqs. fIBTj) and fl68l) self-consistently, and then 
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we have [93 



M n m 



1 



M 



exp [-(3E(x m (k))] 



m=l k=l 



Y,niexp(ft)WW(E(x m (k))) 



< A> T 



M rim 



1 



(69) 



EE 



M 



exp [-pE(x m (k))] 



m=l fc=l 



E"<opOi)^W*))) 



where the trajectories x m (k) (fc = 1, ■ • ■ taken from each multicanonical simula- 

tion with the multicanonical weight factor W^(E) (m — 1, • ■ ■ , A4) separately. 



As seen above, both REMUCA and MUCAREM can be used to obtain the multi- 
canonical weight factor, or the density of states, for the entire potential energy range of 
interest. For complex systems, however, a single REMUCA or MUCAREM simulation is 
often insufficient. In such cases we can iterate MUCA (in REMUCA) and/or MUCAREM 
simulations in which the estimate of the multicanonical weight factor is updated by the 
single- and/or multiple-histogram reweighting techniques, respectively. 

To be more specific, this iterative process can be summarized as follows. The RE- 
MUCA production run corresponds to a MUCA simulation with the weight factor W mu (E). 
The new estimate of the density of states can be obtained by the single-histogram reweight- 
ing techniques of Eq. (]27[) . On the other hand, from the MUCAREM production run, 
the improved density of states can be obtained by the multiple-histogram reweighting 
techniques of Eqs. (fBTjl and (|68|) . 

The improved density of states thus obtained leads to a new multicanonical weight 
factor (see Eq. (1171) ). The next iteration can be either a MUCA production run (as in 
REMUCA) or MUCAREM production run. The results of this production run may yield 
an optimal multicanonical weight factor that yields a sufficiently flat energy distribution 
for the entire energy range of interest. If not, we can repeat the above process by obtaining 
the third estimate of the multicanonical weight factor either by a MUCA production run 
(as in REMUCA) or by a MUCAREM production run, and so on. 

We remark that as the estimate of the multicanonical weight factor becomes more 
accurate, one is required to have a less number of replicas for a successful MUCAREM 
simulation, because each replica will have a flat energy distribution for a wider energy 
range. Hence, for a large, complex system, it is often more efficient to first try MU- 
CAREM and iteratively reduce the number of replicas so that eventually one needs only 
one or a few replicas (instead of trying REMUCA directly from the beginning and iterat- 
ing MUCA simulations). 

We now describe the simulated tempering replica-exchange method (STREM) |115j . 
Suppose that the simulated tempering weight factor Wst(E; T n ) (or equivalently, the di- 
mensionless Helmholtz free energy a n in Eq. has been obtained as in REST or by 

any other methods in the entire temperature range of interest (7\ < T n < Tm)- We divide 
the overlapping temperature ranges into M. regions {M. <C M). Suppose each tempera- 
ture range m has N m temperatures: (fc = 1, • • • , J\f m ) for m — 1, • ■ • , M.. We assign 
each temperature range to a replica; each replica i is in one-to-one correspondence with 



are 
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a different temperature range m of ST run, where T x < < Tjj 1 ' {k — 1, ■ ■ • ,Af m ). 
We then introduce the replica-exchange process between neighboring temperature ranges. 
This works when we allow sufficient overlaps between the temperature regions. 

A STREM simulation is then realized by alternately performing the following two 
steps. |115] 

1. Each replica performs a ST simulation within the fixed temperature range simultaneously 
and independently for a certain MC or MD steps. 

2. A pair of replicas, say i and j, which are at, say T = T^ m ^ and T = T^ m+1 \ in neigh- 
boring temperature ranges, say m-th and (m + l)-th, respectively, are exchanged: 
X = j- • • , x§ , ■ ■ • , xj/', ■ ■ ■} — > X' = j- ■ • , ■ ■ ■ , x$ , ■ ■ -}. The transition prob- 
ability of this replica exchange is given by the Metropolis criterion: 

w(X =min(l,exp(-A)) , (70) 

where 

A=(^-pt^)(E( q M)-E(a®)) . (71) 

While in MUCAREM each replica performs a random walk in multicanonical ensemble 
of finite energy range, in STREM each replica performs a random walk by simulated 
tempering of finite temperature range. These "local" random walks are made "global" to 
cover the entire energy range of interest by the replica-exchange process. 

2.5 Multidimensional Replica- Exchange Method 

We now present our multidimensional extension of REM, which we refer to as multidi- 
mensional replica- exchange method (MREM) [86]. The crucial observation that led to the 
new algorithm is: As long as we have M non-interacting replicas of the original system, 
the Hamiltonian H(q,p) of the system does not have to be identical among the replicas 
and it can depend on a parameter with different parameter values for different replicas. 
Namely, we can write the Hamiltonian for the i-th replica at temperature T m as 

H m (g w ,P w ) = A'(p w ) + ^A m (g w ) , (72) 

where the potential energy E\ m depends on a parameter A m and can be written as 

E Xm (q®) = E (q®)+\ m V(q®) . (73) 

This expression for the potential energy is often used in simulations. For instance, in 
umbrella sampling [21], Eo(q) and V(q) can be respectively taken as the original potential 
energy and the "biasing" potential energy with the coupling parameter A m . In simulations 
of spin systems, on the other hand, E (q) and V(q) (here, q stands for spins) can be 
respectively considered as the zero-field term and the magnetization term coupled with 
the external field A m . 

While replica i and temperature T m are in one-to-one correspondence in the original 
REM, replica i and "parameter set" A m = (T m , A m ) are in one-to-one correspondence 
in the new algorithm. Hence, the present algorithm can be considered as a multidimen- 
sional extension of the original replica-exchange method where the "parameter space" is 
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one-dimensional (i.e., A m = T m ). Because the replicas are non- interacting, the weight 
factor for the state X in this new generalized ensemble is again given by the product of 
Boltzmann factors for each replica (see Eq. (j42p ): 

r M 

Wmrem(X) = exp I -J2 Pm(i)H m{i) (q® , p w ) 

= exp {-£ P m H m (gt*M],p[<(«)]) 1 , 
I m=l J 

where i(m) and m(i) are the permutation functions in Eq. (T4T31) . Then the same derivation 
that led to the original replica-exchange criterion follows, and the transition probability 
of replica exchange is given by Eq. (j551) . where we now have (see Eq. floUl) ) [86] 

A = An (g b] ) - ^ m (g w )) - A. K (g b1 ) - ^ (g [i] )) • (75) 

Here, i?A m and E\ n are the total potential energies (see Eq. (1731) ). Note that we need to 
newly evaluate the potential energy for exchanged coordinates, E\ m (q^) and E\ n (q^), 
because E\ m and E\ n are in general different functions. 

For obtaining the canonical distributions, the multiple-histogram reweighting tech- 
niques [HI [HI] are particularly suitable. Suppose we have made a single run of the present 
replica-exchange simulation with M replicas that correspond to M different parame- 
ter sets A m = (T m ,X m ) (m = 1, ■■•,M). Let N m (E , V) and n m be respectively the 
potential-energy histogram and the total number of samples obtained for the m-th pa- 
rameter set A m . The WHAM equations that yield the canonical probability distribution 
Pt,\(Eq,V) = n(E ,V) exp(—j3E\) with any potential-energy parameter value A at any 
temperature T = 1/ksP are then given by [86] 

M 

J2N m {E ,V) 

n(E , V) = — m=1 , (76) 

n m exp (f m - (3 m E Xrn ) 

m=l 

and for each m (= 1, • ■ ■ , M) 

exp(-/ m ) = 5] n{E , V) exp {-(3 m E Xm ) . (77) 

Eo,V 

Here, n(Eo,V) is the generalized density of states. Note that n(E Q ,V) is independent 
of the parameter sets A m = (T m , A m ) (m = 1,---,M). The density of states n(E ,V) 
and the "dimensionless" Helmholtz free energy f m in Eqs. ( !76|) and (ITT)) are solved self- 
consistently by iteration. 

We can use MREM for free energy calculations. We first describe the free-energy 
perturbation case. The potential energy is given by 

EM = E^q) + A (E F (q) - EM)) , (78) 

where Ej and E F are the potential energy for a "wild-type" molecule and a "mutated" 
molecule, respectively. Note that this equation has the same form as Eq. (I75|) . 
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Our replica-exchange simulation is performed for M replicas with M different values 
of the parameters A m = (T m ,A m ). Since E\ =0 (q) = Ej(q) and E\=i(q) = E F (q), we 
should choose enough A m values distributed in the range between and 1 so that we 
may have sufficient acceptance of replica exchange. From the simulation, M histograms 
N m (Ej,EF — Ej), or equivalently N m (Ej, Ep), are obtained. The Helmholtz free energy 
difference of "mutation" at temperature T (= 1//cb/3), AF = F\=i — F\ = q, can then be 
calculated from 

£ P T ^ =l {E u E F ) 
ex P (-/3AF) = Z f^ = — - , (79) 

Ej,Ep 

where Px t x(Ej, Ep) = n(Ej, Ep) exp (—f3E\) are obtained from the WHAM equations of 
Eqs. (j75j) and (TTTj). 

We now describe another free energy calculations based on MREM applied to umbrella 
sampling [21], which we refer to as replica- exchange umbrella sampling (REUS). The 
potential energy is a generalization of Eq. ( 1731) and is given by 

E\(q)=E (Q)+ y t*®V t (q) , (80) 

where E (q) is the original unbiased potential, Vi(q) (£ — 1, ■ ■ ■ , L) are the biasing (um- 
brella) potentials, and are the corresponding coupling constants (A = (\^\ ■ • • , A^)). 
Introducing a "reaction coordinate" £, the umbrella potentials are usually written as har- 
monic restraints: 

Vt(q)=k e (t(q)-d t ) 2 , (£ = 1,---,L) , (81) 

where di are the midpoints and kg are the strengths of the restraining potentials. We 
prepare M replicas with M different values of the parameters A m = (T m ,X m ), and the 
replica-exchange simulation is performed. Since the umbrella potentials Ve(q) in Eq. ( 1HT1) 
are all functions of the reaction coordinate £ only, we can take the histogram N m (Eo,£) 
instead of N m (E , Vi, ■ ■ ■ , Vl). The WHAM equations of Eqs. (1751) and (177)) can then be 
written as [86] 

M 

E N rn (E ,0 

n(Eo, = M m=1 (82) 



M 

E rim exp (/ m - f3 m E^ 



m=l 



and for each m (= 1, • ■ ■ , M) 



exp(-/ m ) = n ( E o, exp (-PrnE Xm ) . (83) 

The expectation value of a physical quantity A with any potential-energy parameter value 
A at any temperature T (= l/ksP) is now given by 

$>(£ ,OP TiA (£o,0 
<A> rA = ^- , (84) 
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where P t ^(Eq,£) = n(E , £) exp (— j3E^j is obtained from the WHAM equations of 
Eqs. (JH21) and (|83|). 

The potential of mean force (PMF), or free energy as a function of the reaction coor- 
dinate, of the original, unbiased system at temperature T is given by 



Wt,A={o}(0 = -^bT1u 



55) 



where {0} = (0, •••,0). 

We now present two examples of realization of REUS. In the first example, we use only 
one temperature, T, and L umbrella potentials. We prepare replicas so that the potential 
energy for each replica includes exactly one umbrella potential (here, we have M = L). 
Namely, in Eq. (IHUi) for A = A m we set 

A$ = km , (86) 

where 5k,i is Kronecker's delta function, and we have 

E Kn (q®) = E {q®) + V m {q®) . (87) 

We exchange replicas corresponding to "neighboring" umbrella potentials, V m and V m+ i. 
The acceptance criterion for replica exchange is given by Eq. (1531) . where Eq. (!75|) now 
reads (with the fixed inverse temperature j3 = l/k^T) [86] 

A = P (V m (qW) - V m (gW) - V m+1 (qW) + V m+1 (f)) , (88) 

where replicas % and j respectively have umbrella potentials V m and V m+ i before the 
exchange. 

In the second example, we prepare Nt temperatures and L umbrella potentials, which 
makes the total number of replicas M = Nt x L. We can introduce the following re- 
labeling for the parameters that characterize the replicas: 

A m = (T m , A m ) — > A/ 5 j = (Tf , Aj) . , s 

(m = l,---,M) (/=!,-••, N T , J=1,---,L) { ™> 



The potential energy is given by Eq. flHTj) with the replacement: m — > J. We perform the 
following replica-exchange processes alternately: 

1. Exchange pairs of replicas corresponding to neighboring temperatures, Tj and X/ + i 
(i.e., exchange replicas i and j that respectively correspond to parameters A/ j and 
A 

i+i,j)- (We refer to this process as T-exchange.) 

2. Exchange pairs of replicas corresponding to "neighboring" umbrella potentials, Vj 
and Vj+i (i.e., exchange replicas % and j that respectively correspond to parameters 
A/ t j and A^j+i). (We refer to this process as A-exchange.) 

The acceptance criterion for these replica exchanges is given by Eq. (155]) . where Eq. (1751) 
now reads [86] 

A = (J3j - (3 I+1 ) (Eo (q b] ) + Vj (qM) - Eo (q [t] ) - Vj (f)) , (90) 

for T-exchange, and 

A = (Vj (gM) - Vj (gW) - V J+1 (gW) + V J+1 (gW)) , (91) 

for A-exchange. By this procedure, the random walk in the reaction coordinate space as 
well as in the temperature space can be realized. 
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2.6 From Multidimensional REM to Multidimensional MUCA 
and ST 

The formulations of MREM give multidimensional/multivariable extensions of REMUCA 
and REST [5]. In REMUCA and in REST, the multicanonical weight factor and the simu- 
lated tempering weight factor are determined from the results of a short REM simulation, 
respectively. The results of a short MREM simulation can therefore be used to determine 
the weight factors for multidimensional/multivariable MUCA and ST simulatoins, where 
random walks in multidimensional "energy" and "parameter" space are realized [5]. Here, 
we give more details. 

We consider a simple example with the following potential energy: 

E x (q) = E (q) + XV(q) . (92) 

In the two-dimensional multicanonical ensemble each state is weighted by the multicanon- 
ical weight factor W mu (E , V) so that a uniform potential energy distribution both in E 
and V may be obtained: 

P nm (Eo, V) oc n{E , V)W mu (E , V) = const , (93) 

where n(E , V) is the two-dimensional density of states. This implies that 

W mn (E , V) = exp [-/3 £ m u(£o, V; T )} = 1 , (94) 

n{& , v ) 

where we have chosen an arbitrary reference temperature, T = l/k^fio, and the "multi- 
canonical "potential energy" is defined by 

E mn {E , V- T ) ee k B T \nn(E , V) . (95) 

The two-dimensional MUCA MC simulation can be performed with the following tran- 
sition probability from state x with potential energy E + XV to state x' with potential 
energy E ' + XV (see Eq. (fT9|)): 

w(x — ► x ) — mm 1, — — — — — — — = mm 1, ——. T r , N . 96 

The MD algorithm in the two-dimensional multicanonical ensemble also naturally follows 
from Eq. ( |T7|) . in which the regular constant temperature MD simulation (with T = To) 
is performed by replacing E by E mu in Eq. ( |T2l) (see Eq. ( l2Til ): 

dE mu (E ,V;T ) s 
Pk = W k ~s Pk ' (97) 

In the two-dimensional simulated tempering, the parameter set (T, A) become dynam- 
ical variables, and both the configuration and the parameter set are updated during the 
simulation with a weight (see Eq. (]30p ): 

W ST (E X] T, X) = exp (—/3E\ + f(T, A)) , (98) 
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where the function /(T, A) is chosen so that the probability distribution of the two- 
dimensional parameter set is flat (see Eq. (1311) ): 

P ST (T, A) = / dE dV n(E , V) W ST (E x] T, A) 

= / dE dV n(E 0l V) exp (-/3E X + f(T, A)) = const . { } 

In the numerical work we discretize the parameter set in M = Nt x L different values, 
(T/ , Aj) (/ = 1, • • • , Nt, J = 1, ■■•,£). Without loss of generality we can order the 
parameters so that T\ < T 2 < • • • < Tn t and Ai < A 2 < • • ■ < Xl- The free energy / is 
now written as fi t j = f(Tj, Aj). Once the initial configuration and the initial parameter 
set are chosen, the two-dimensional ST is then realized by alternately performing the 
following two steps: 

1. A canonical MC or MD simulation at the fixed parameter set (Xj, Aj) is carried out 
for a certain steps. 

2. One of the parameters in the parameter set (Tj, Aj) is updated to the neighbor- 
ing values with the configuration and the other parameter fixed. The transition 
probability of this parameter-updating process is given by the following Metropolis 
criterion: 

w(Tj - T I±1 ) = min ( 1, ^f"^^ ) = min (1, exp (-A)) , (100) 

\ w ST {hxj; J-i, Aj) / 



where 

for T-update, and 



A = ((3 I±1 - Pj) E Xj - (f I±ljJ - f ItJ ) , (101) 



w{Xj -> A J±1 ) = mm 1, — - = mm (1, exp (-A)) , (102) 



where 



for A-update. 



A =Pj (E Xj±1 - E Xj ) - (/ J)J±1 - f^j) 
= Pi{\j±i - h)V - (fi,J±i - fi,j) 



(103) 



Finally, we present the corresponding MREM. We prepare Nt temperatures and L 
A parameters, which makes the total number of replicas M = Nt x L. We perform the 
following replica-exchange processes alternately: 

1. Exchange pairs of replicas corresponding to neighboring temperatures, Tj and Tj+i 
(We refer to this process as T-exchange.) 

2. Exchange pairs of replicas corresponding to "neighboring" A parameters, Aj and 
A j+i (We refer to this process as A-exchange.) 

The acceptance criterion for these replica exchanges is given by Eq. (153j) . where Eq. ( 1751) 
now reads 

A = ((3j - (3 I+1 ) (E Xj ( g W) - E Xj ( g W)) , (104) 
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for T-exchange, and 

A = 0j (E Xj (gW) - E Xj (gW) - £ Aj+1 (gW) + E Xj+1 (q 

= A(A J -A J+1 )(y(gM)-y(gW)) , (105) 

for A-exchange. 

After a short MREM simulation, we can use the multiple-histogram reweighting tech- 
niques to obtain n(E , V) and fij. Let Ni t j(E , V) and nij be respectively the potential- 
energy histogram and the total number of samples obtained for the parameter set (Tj, Xj). 
The WHAM equations are then given by 

N T L 

E E ^(£ , v) 
E E n i,j ex p (.fifj - PiExj) 

1=1 J=l 

and for each / and J (I — 1, • • ■ , N T , J = 1, ■ ■ • , L) 

exp (-fjj) = Y, n{E , V) exp {-^E Xj ) . (107) 

Eq,V 

These equations are solved self-consistently by iteration for n(E , V) and fi t j. 

Hence, we can determine the multidimensional multicanonical weight factor W mn (Eo, V) 
and the multidimensional simulated tempering weight factor Wst(E\/, Ti, Xj). The for- 
mer is given by 

W mu (E ,V)= 1 , (108) 
n(E , V) 

and the latter is given by 



W ST (E Xj ; Tj, Xj) = exp (-PjE Xj + f ItJ ) . (109) 



3 SIMULATION RESULTS 

We first compare the performances of REM, MUCAREM, and REMUCA. The accuracy of 
average quantities calculated depend on the "quality" of the random walk in the potential 
energy space, and the measure for this quality can be given by the number of tunneling 
events [151 US]- O ne tunneling event is defined by a trajectory that goes from En to E^ 
and back, where En and E^ are the values near the highest energy and the lowest energy, 
respectively, which the random walk can reach. If En is sufficiently high, the trajectory 
gets completely uncorrelated when it reaches En- On the other hand, when the trajectory 
reaches near E^, it tends to get trapped in local-minimum states. We thus consider that 
the more tunneling events we observe during a fixed number of MC/MD steps, the more 
efficient the method is as a generalized-ensemble algorithm (or, the average quantities 
obtained by the reweighting techniques are more reliable). 

The first example is Monte Carlo simulations of the system of a 17-residue fragment 
of ribonuclease Tl in implicit solvent (expressed by the solvent accessible surface area) 
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[94] . The amino-acid sequence is Ser-Ser-Asp-Val-Ser-Thr-Ala-Gln-ile-Ala-Ala-Tyr-Lys- 
Leu-His-Glu-Asp. The energy function E^ot that we used is the sum of the confor- 
mational energy term of the solute Ep and the solvation free energy term E^ol for the 
interaction of the peptide with the surrounding solvent: -Etot = Ep + -Esol- Here, 
the solvation term -Esol is given by the sum of the terms that are proportional to the 
solvent-accessible surface area of the atomic groups of the solute. The parameters in 
the conformational energy as well as the molecular geometry were taken from ECEPP/2. 
The parameters of the solvent term were adopted from Ref. [125j . The computer code 
KONF90 [1261 H23 was used, and MC simulations based on the REM, MUCAREM, and 
REMUCA were performed. For the calculation of a solvent-accessible surface area, we 
used the computer code NSOL |128j . The dihedral angles <fi and ip in the main chain 
and x i n the side chain constituted the variables to be updated in the MC simulations. 
The number of degrees of freedom for the peptide is 80. One MC sweep consists of up- 
dating all these angles once with Metropolis evaluation for each update. The simulations 
were started from randomly generated conformations. In Table [1] we list the number of 
tunneling events in REM, MUCAREM, and REMUCA simulations of the same system 

m- 

Table 1: Number of tunneling events in the MC simulations of a fragment of ribonuclease 
Tl for REM, MUCAREM, and REMUCA simulations 



Total MC sweeps 


REM 


MUCAREM 


REMUCA 


2 x 10 6 


2 


9 


18 


3 x 10 6 


5 


16 


29 


4 x 10 6 


9 


22 


38 



Table 2: Number of tunneling events in the MD simulations of three peptides in expicit 
water for REM, MUCAREM, and REMUCA simulations 



Peptide 


No. of atoms 


Total MD steps 


REMD 


MUCAREM 


REMUCA 


Alanine dipeptide 


418 


4 x 10 6 


11 


40 


59 


Alanine trimer 


876 


5 x 10 6 


1 


20 


29 


Met-enkephalin 


1662 


8 x 10 6 





12 


27 



Hence, REMUCA is the most efficient, then MUCAREM, and finally REM. 

The next systems are small peptides in explicit water |129j . When we consider explicit 
water molecules, the problem becomes order-of-magnitude more difficult than the case 
with implicit water models. They are alanine dipeptide with 132 water molecules, alanine 
trimer with 278 water molecules, and Met-enkephalin with 526 water molecules. The 
force-field, or the potential energy, that we used is AMBER parm96 [130] for the peptides 
and TIP3P |131j for water molecules. The peptides were placed inside the spheres of 
water molecules and the harmonic constraining forces were imposed in order to prevent 
the water molecules from flying apart. The unit time step, At, was set to 0.5 fsec. The 
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modified version [1321 1133] of the software PRESTO version 2 [134] was used. In Table [2] 
we list the number of tunneling events in these systems. 

The last system is the C-peptide of ribonuclease A in explicit water [135j . In the 
model of simulations, the N-terminus and the C-terminus of the C-peptide analogue were 
blocked with the acetyl group and the N-methyl group, respectively. The number of 
amino acids is 13 and the amino-acid sequence is: Ace-Ala-Glu~-Thr-Ala-Ala-Ala-Lys + - 
Phe-Leu-Arg + -Ala-His + -Ala-Nme [1361 1137] . The initial configuration of our simulation 
was first generated by a high temperature molecular dynamics simulation (at T = 1000 
K) in gas phase, starting from a fully extended conformation. We randomly selected one 
of the structures that do not have any secondary structures such as a- helix and /3-sheet. 
The peptide was then solvated in a sphere of radius 22 A, in which 1387 water molecules 
were included (see Fig. 1). Harmonic restraint was applied to prevent the water molecules 
from going out of the sphere. The total number of atoms is 4365. The dielectric constant 
was set equal to 1.0. The force-field parameters for protein were taken from the all-atom 
version of AMBER parm99 [138] , which was found to be suitable for studying helical 
peptides |139] . and TIP3P model [13 lj was used for water molecules. The unit time step, 
At, was set to 0.5 fsec. 

In Table [3] the essential parameters in the simulations performed in this article are 
summarized. 

We first performed a REMD simulation with 32 replicas for 100 psec per replica 
(REMD1 in Table [3]). During this REMD simulation, replica exchange was tried every 
200 MD steps. Using the obtained potential-energy histogram of each replica as input data 
to the multiple-histogram analysis in Eqs. (4) and (5), we obtained the first estimate of 
the multicanonical weight factor, or the density of states. We divided this multicanonical 
weight factor into four multicanonical weight factors that cover different energy regions 
[88| l93l 194"] and assigned these multicanonical weight factors into four replicas (the weight 
factors cover the potential energy ranges from —13791.5 to —11900.5 kcal/mol, from 
-12962.5 to -10796.5 kcal/mol, from -11900.5 to -9524.5 kcal/mol, and from -10796.5 
to —8293.5 kcal/mol). We then carried out a MUCAREM simulation with four replicas 
for 1 nsec per replica (MUCAREM1 in Table [3]), in which replica exchange was tried every 
1000 MD steps. We again used the potential-energy histogram of each replica as the input 
data to the multiple-histogram analysis and finally obtained the multicanonical weight 
factor with high precision. As a production run, we carried out a 15 nsec multicanonical 
MD simulation with one replica (REMUCAl in Table [3D and the results of this production 
run were analyzed in detail. 

In Fig. [2] we show the probability distributions of potential energy that were obtained 
from the above three generalized-ensemble simulations, namely, REMD1, MUCAREM1, 
and REMUCAl. We see in Fig. [2](a) that there are enough overlaps between all pairs 
of neighboring canonical distributions, suggesting that there were sufficient numbers of 
replica exchange in REMD1. We see in Fig. W(Jo) that there are good overlaps between 
all pairs of neighboring multicanonical distributions, implying that MUCAREM1 also 
performed properly. Finally, the multicanonical distribution in Fig. [2](c) is completely flat 
between around —13000 kcal/mol and around —8000 kcal/mol. The results suggest that 
a free random walk was realized in this energy range. 

In Fig. [3h we show the time series of potential energy from REMUCAl. We indeed 
observe a random walk covering as much as 5000 kcal/mol of energy range (note that 23 
kcal/mol ~ 1 eV). We show in Fig. [3^b) the average potential energy as a function of 



26 



temperature, which was obtained from the trajectory of REMUCAl by the reweighting 
techniques. The average potential energy monotonically increases as the temperature 
increases. 

Here, we took Eh = —8250 kcal/mol and = —12850 kcal/mol for the measurement 
of the tunneling events. The random walk in REMUCAl yielded as many as 55 tunneling 
events in 15 nsec. The corresponding numbers of tunneling events for REMD1 and for 
MUCAREM1 were in 3.2 nsec and 5 in 4 nsec, respectively. Hence, REMUCA is the 
most efficient and reliable among the three generalized-ensemble algorithms. 

In Fig.H]the potential of mean force (PMF), or free energy, along the first two principal 
component axes at 300 K is shown. There exist three distinct minima in the free-energy 
landscape, which correspond to three local-minimum-energy states. We show represen- 
tative conformations at these minima in Fig. [51 The structure of the global-minimum 
free-energy state (GM) has a partially distorted a-helix with the salt bridge between 
Glu~-2 and Arg + -10. The structure is in good agreement with the experimental struc- 
ture obtained by both NMR and X-ray experiments. In this structure there also exists a 
contact between Phe-8 and His + -12. This contact is again observed in the corresponding 
residues of the X-ray structure. At LM1 the structure has a contact between Phe-8 and 
His + -12, but the salt bridge between Glu~-2 and Arg + -10 is not formed. On the other 
hand, the structure at LM2 has this salt bridge, but it does not have a contact between 
Phe-8 and His + -12. Thus, only the structures at GM satisfy all of the interactions that 
have been observed by the X-ray and other experimental studies. 

Finally, we remark that the largest peptide in explicit water that we have succeeded in 
folding into the native structure from random initial conformations is so far the 16-residue 
C-terminal /3-hairpin of streptococcal protein G Bl domain, which was accomplished by 
MUCAREM simulations with eight replicas [142] . 

4 CONCLUSIONS 

In this article we have reviewed some of powerful generalized-ensemble algorithms for 
both Monte Carlo simulations and molecular dynamics simulations. A simulation in 
generalized ensemble realizes a random walk in potential energy space, alleviating the 
multiple-minima problem that is a common difficulty in simulations of complex systems 
with many degrees of freedom. 

Detailed formulations of the three well-known generalized-ensemble algorithms, namely, 
multicanonical algorithm (MUCA), simulated tempering (ST), and replica-exchange method 
(REM), were given. 

We then introduced several new generalized-ensemble algorithms that combine the 
merits of the above three methods. 

The question is then which method is the most recommended. Our criterion for 
the effectiveness of generalized-ensemble algorithms was how many random walk cycles 
(tunneling events) in potential energy space between the high-energy region and low- 
energy region are realized within a fixed number of total MC (or MD) steps. We found 
that once the optimal MUCA weight factor is obtained, MUCA (and REMUCA) is the 
most effective (i.e., has the most number of tunneling events), and REM is the least 
[93] . We also found that once the optimal ST weight factor is obtained, ST (and REST) 
has more tunneling events than REM [89l I115j . Moreover, we compared the efficiency 
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of Berg's recursion [73], Wang-Landau method [231 [21], and REMUCA/MUCAREM as 
methods for the multicanonical weight factor determination in two-dimensional 10-state 
Potts model and found that the three methods are about equal in efficiency |143j - [T4"5] . 

Hence, the answer to the above question will depend on how much time one is willing 
to (or forced to) spend in order to determine the MUCA or ST weight factors. Given a 
problem, the first choice is REM because of its simplicity (no weight factor determination 
is required). If REM turns out to be insufficient or too much time-consuming (like the 
case with first-order phase transitions), then other more powerful algorithms such as 
MUCAREM and STREM are recommended. 
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Figure 1: The initial configuration of C-peptide in explicit water, which was used in all of 
the 32 replicas of the first REMD simulation (REMD1 in Tabled]). The red filled circles 
stand for the oxygen atoms of water molecules. The number of water molecules is 1387, 
and they are placed in a sphere of radius 22 A. As for the peptide, besides the backbone 
structure (in blue), side chains of only Glu~-2, Phe-8, Arg + -10, and His + -12 are shown 
(in yellow). The figure was created with Molscript [140] and Raster3D |141j . 



Table 3: Summary of parameters in REMD, MUCAREM, and REMUCA simulations 





Number of 


Temperature, 






MD steps per 




replicas, M 


T m (K) (m 


= 1, 


••,M) 




replica 


REMD1* 


32 


250, 


258, 267, 


276, 


286, 295, 


305, 


2.0 x 10 5 






315, 


326, 337, 


348, 


360, 372, 


385, 








398, 


411, 425, 


440, 


455, 470, 


486, 








502, 


519, 537, 


555, 


574, 593, 


613, 








634, 


655, 677, 


700 








MUCAREM1 


4 


360, 


440, 555, 


700 






2.0 x 10 6 


REMUCA1 


1 


700 










3.0 x 10 7 



* REMD1 stands for the replica-exchange molecular dynamics simulation, MUCAREM1 
stands for the multicanonical replica-exchange molecular dynamics simulation, and RE- 
MUCAl stands for the final multicanonical molecular dynamics simulation (the produc- 
tion run) of REMUCA. The results of REMD1 were used to determine the multicanonical 
weight factors for MUCAREM1, and those of MUCAREM1 were used to determine the 
multicanonical weight factor for REMUCAl. 
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Figure 2: Probability distributions of potential energy of the C-peptide system obtained 
from (a) REMD1, (b) MUCAREM1, and (c) REMUCAl. See Table [3]for the parameters 
of the simulations. Dashed curves in (c) are the reweighted canonical distributions at 290, 
300, 500, and 700 K (from left to right). 
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Figure 3: Time series of potential energy of the C-peptide system from the REMUCA 
production run (REMUCAl in Table [3]) (a) and the average potential energy as a function 
of temperature (b). The latter was obtained from the trajectory of REMUCAl by the 
single-histogram reweighting techniques. 
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Figure 4: Potential of mean force (kcal/mol) of the C-peptide system along the first 
two principal components at 300 K. The free energy was calculated from the results of 
REMUCA production run (REMUCAl in Table [3]) by the single-histogram reweighting 
techniques and normalized so that the global-minimum state (GM) has the value zero. 
GM, LM1, and LM2 represent three distinct minimum free-energy states. 



a b □ 

Figure 5: The representative structures at the global-minimum free-energy state ((a) GM) 
and the two local- minimum states ((b) LM1 and (c) LM2). As for the peptide structures, 
besides the backbone structure, side chains of only Glu~-2, Phe-8, Arg + -10, and His + -12 
are shown in ball-and-stick model. 
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